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1 Introduction 

With the steadily increasing computing power multivariate methods are nowadays standard tech¬ 
niques in the experimental analysis. Initially introduced in experimental studies where the event 
rates are small and signals are difficult to disentangle from overwhelming backgrounds, multivariate 
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methods have been proven useful also in a more general context. Among these methods the Matrix 
Element Method represents a prominent example, since it allows a direct comparison of observed 
event samples with expectations within a specific theoretical model. As originally introduced in 
refs. [1,2] the method is based on the assumption that the probability to observe a specific event 
can be calculated using the corresponding matrix element together with so-called transfer functions 
which model the probability for the observation of a given partonic event as a specific (hadronic) 
event at the detector level. It has been argued that if all the ingredients in this procedure are known 
with optimal accuracy, the event likelihood defined in this way represents an optimal statistical test. 
The method thus makes maximal use of the information contained in the single event. Based on 
this assumption the Matrix Element Method has been used for example to measure the top-quark 
mass at the Tevatron (see for example refs. [3-5]). Alternatively the Matrix Element Method can 
be used to distinguish different hypotheses like for example background versus signal hypothesis 
or SM versus BSM physics hypothesis for a given event sample (see for example ref. [6, 7]). Eet 
us also mention that the Matrix Element Method exists today in different flavors like for example 
the MEEA approach as used for example in refs. [8, 9]. The construction of optimal observables 
as used for example in refs. [10-12] may also be seen as a variation of the same theme. So far in 
most cases the method is only applied using leading-order (EO) matrix elements. To simplify the 
application of the Matrix Element Method the automated calculation of the required event weights 
has been studied recently in ref. [13]. Given the recent progress concerning the calculation of NEO 
corrections it is natural to include also NEO corrections in the evaluation of the matrix elements. A 
first attempt in this direction has been made in ref. [14] where the effect of QCD radiation has been 
studied. In refs. [15, 16] the radiation pattern of boosted Higgs bosons and top quarks is studied 
and compared with the radiation profile of QCD jefs. In ref. [17] the information from the hard 
matrix element and a parton shower is used for a signal versus background discrimination for the 
signal process Z' decaying to boosted top quarks. In refs. [18-20] the impact of NEO corrections 
including also the virtual corrections is investigated and a possible extension of the Matrix Ele¬ 
ment Method beyond the Born approximation is discussed. A first detailed application has been 
presented in ref. [21] where Higgs production with subsequent decay into H ^ Zy is investigated. 
So far the method presented in refs. [18, 19] is restricted to the production of uncolored objects 
and the extension to the production of colored particles like for example top-quark pairs is still 
missing. In ref. [20] the extension to include hadronic production of jets is investigated by means 
of a longitudinal boost along the beam axis to remove the unbalanced transverse momentum and 
map NEO and EO jets. 

The extension of the final slate phase space encountered in the generalisation of the Matrix Element 
Method beyond the Bom approximation is an intrinsic problem which makes the NEO extension 
non-trivial. Real corrections which appear as additional contribution, when higher order correc¬ 
tions are taken into account, allow for the emission of an additional parton. In addition to the 2^n 
Bom kinematics also contributions living on an n -i-1 parton phase space need thus to be considered. 
Restricting the attention to jet physics, one may argue that only regions of the extended phase space 
in which the additional parton is clustered/merged into a jet contribute to n-jet observables and we 
end up again with a 2 ^ n configuration—now in terms of jets instead of partons. However, ap¬ 
plying standard jet algorithms, the resulting jets typically do not satisfy the kinematical constraints 
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of the Bom process. In particular, the clustering will in general create a non vanishing mass for 
the jets even for massless partons. Also momentum conservation is not necessarily required by all 
jet algorithms. It is thus not clear how the contribution from the real corrections can be unambigu¬ 
ously combined with the virtual corrections which respect the Bom kinematics. Furthermore, the 
practical calculation of next-to-leading order event weights—which must take into account virtual 
as well as real corrections—is non-trivial. As we shall describe in more detail in the next section, 
the problem is related to the phase space integration of real corrections which lead to an n-particle 
final state after clustering. Defining event weights for ‘jet-events’ in NLO accuracy thus requires 
two aspects to be addressed: 

1. In the theoretical predictions a modified clustering may be used which guarantees momentum 
conservation and keeps the clustered jets on-shell. More precisely, the mass of a jet which 
is formed through a merging of a massless parton with another parton should be equal to 
the mass of the radiating parton. Merging two massless partons should result in a massless 
jet. It is worth stressing that also for established jet algorithms, which rely on a 2 —> 1 
clustering, the jet masses created in the perturbative calculation have very little to do (in 
case of massless partons) with the jet masses observed in the experimentally measured jets 
in case of ‘light’ jets.^ The latter are mostly related to non-perturbative effects. In that sense, 
a modified clustering as proposed here is equally well motivated as what is currently used. 
In fact, one may even argue that a clustering fulfilling the above constraints may lead to a 
better separation of perturbative and non-perturbative physics. 

2. A method needs to be constructed allowing the efficient integration of the regions in the 
n + 1-parton phase space contributing to the considered n-jet configuration. (In principle an 
inclusive observable may also receive contributions from n -i-1 jet configurations requiring 
the evaluation of the corresponding weights. In this case the problem is however very similar 
to the leading-order situation.) 

In this article we illustrate that both aspects can be addressed by using an appropriate jet clustering 
which is intimately related to a factorised parameterisation of the n + 1-parton phase space into an 
n-particle phase space of the n jets times the phase space related to the ‘clustered’ (unobserved) 
parton. The general idea is to extend the typical 2 —> 1 clustering to a 3 —> 2 clustering as it is 
well known for example from the dipole subtraction method [22, 23]. This clustering satisfies 
momentum conservation as well as on-shell conditions at the expense of introducing an additional 
spectator which allows to guarantee momentum conservation which would be otherwise violated 
by enforcing the on-shell condition. Some freedom exists how to choose the additional ‘spectator’. 
For example, to minimise the difference to the traditional clustering, one may choose the spectator 
such that the momentum reshuffling is minimised. Having chosen the spectator, a recombination 
according to the Catani-Seymour phase space factorisation [22, 23] is applied. In this way, the 
reduced kinematics appearing in the Catani-Seymour subtraction formalism can be identified with 
the final state jets. At the same time, the factorised phase space can be integrated over the unre- 

*This statement does not apply to fat jets or highly boosted objects. 
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solved regions to obtain the contributions from the real corrections to the event weight. We will 
give more details in the next section. 

It should be noted that using the four mappings given in the Catani-Seymour subtraction algorithm 
as clustering prescriptions in the proposed 3 —> 2 jet algorithm is sufficient, to construct appropriate 
jet algorithms covering most of the relevant collider experiments: production of at least 2 jets at a 
lepton collider, deep inelastic scattering and hadronic production of electroweak final states and/or 
jets. The respective final states can have arbitrary masses. 

Let us briefly compare the method outlined above with some existing work where similar ideas have 
been applied in a different context. For example, the ideas presented here share some features with 
ref. [24]. There are however important differences. In ref. [24] an additional resolution parameter 
is introduced to define ‘resolved’ partons similar to what is done in the phase space slicing method 
[25, 26]. Using the resolved partons in an intermediate step any physical jet algorithm should 
in principle be applicable. Care has to be taken that the two cuts—the artificial one to define 
resolved partons and the physical one to define jets within a given jet algorithm—do not interfere. 
In the approach described above we work directly with jets defined by the jet algorithm used in the 
experimental analysis. In ref. [24] only final state singularities in e^e~ annihilation are considered. 
Here we include initial state singularities as well. Technically, only one phase space mapping is 
required in ref. [24]. As we will see in the next section, in general more mappings are required 
and the application of the Catani-Seymour subtraction method is less obvious. In refs. [27, 28] the 
implementation of a parton shower based on the Catani-Seymour subtraction method is studied. 
Although the aim of this work is rather different, the parameterisation of the phase space is similar 
to what is used in this article and many useful results which we collect in the following can be 
found also in refs. [27, 28]. 

In ref. [29] a method to generate the phase space of n -i-1 massless particles by forward branching 
of configurations of n massless particles is presented. This method which is applied in refs. [18- 
21] employs two 3 —> 2 prescriptions to cluster three massless partons to two massless jets. The 
method that we present in this article can be seen as a generalisation of this approach with two 
further prescriptions for the clustering and the extension to massive particles. As a proof of concept, 
especially for the aforementioned generalisations, we apply the Matrix Element Method in NLO 
accuracy to a process with massive colored particles in the final state. 

The article is organised as follows. In section 3 we present the phase space parameterisation used in 
the numerical integration. Roughly speaking the Catani-Seymour mapping is inverted. In section 

4 we validate the approach by various cross checks. As a proof of concept we illustrate in section 

5 the Matrix Element Method including NEO corrections applied to top-quark pair production in 
e^e~ annihilation. We summarise our main findings in the conclusion in section 6. 
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2 Formalism 


2.1 Matrix Element Method and event weights at next-to-leading order 

The Matrix Element Method tries to make maximal use of the information provided hy an indi¬ 
vidual event. Instead of considering distributions, calculated for event samples, the probability of 
the event in the context of a given theory is investigated. In what follows we assume that all ex¬ 
perimentally available information of the event is collected in the variable x. In the ideal situation 
that all momenta have been reconstructed, one may think of x as the collection of the observed 
momenta which we label with 7i,..., 


(2.1) 

However, since some particles may escape detection or are only partially reconstructed, the exper¬ 
imentally accessible information may be in practice only a subset of this information. 

Putting aside for the moment higher order corrections, one may interpret the partonic cross section 
calculated for a specific model—for example the Standard Model—as the probability distribution to 
observe a partonic event. A model-dependent likelihood, with model parameters Q, for observing 
an event x is than given schematically by 

fdyi...dyn xW{x,f), (2.2) 

crj dyi...dyn 

where the differential cross section is denoted by d"crfdy \... <iy„ and the so-called transfer function 
W{x,y) describes the probability that a partonic event y is measured in the detector as the event x. 
In principle, the variables collected in y may be chosen independently from the variables in x. Even 
the dimension of the two vectors does not need to agree. However, it may prove beneficial to chose 
the two sets as closely related as possible. Assuming that the two can be identified, an ideal detector 
would than correspond to the situation in which the transfer function is given by a delta function: 
Way) = S{^-f). 

Roughly speaking, maximising the likelihood with respect to Q for a given event sample gives an 
estimator for the model parameter. This is the essence of the so-called Matrix Element Method 
(MEM). (More details can be found in refs. [1-3, 5, 7, 13, 30-35].) Since all the information 
available in the single measurement is retained, this approach is believed to make maximal use of 
the information content of the single event. 

While in principle the integration over the transfer functions looks straightforward, in practice it is 
not trivial due to the peak structure of the transfer functions. In addition, we note that the transfer 
functions need to be determined within the experimental analysis, which may also represent a non¬ 
trivial task. In the following we do not consider this issue any further since current experimental 
analyses using the Matrix Element Method are already used to this type of problem. The focus 
of this article is the extension of the MEM beyond the Bom approximation. To be specific we 
assume y in the following to be the collection of final stale momenta—which may be obtained 
in the case of the real corrections through the merging of collinear or soft partons according to a 
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specific clustering procedure. To distinguish these momenta from the partonic momenta we call 
them jet-momenta in what follows. They may be seen as the perturbative approximation of the jets 
observed in the experiments. (More details will be given in the next subsection.) The motivation to 
use the jet-momenta is threefold: 

1. Being differential in the jet-momenta all relevant information about the differential cross 
section is kept, allowing in a second step also to use a different set of variables. 

2. Using the jet momenta, identifying the transfer functions with delta-functions may provide a 
reasonable first approximation. 

3. Being able to calculate event weights for ‘jet events’ including higher order corrections is 
interesting on its own right. 

In the rest of this subsection we illustrate the main obstacle in the calculation of event weights for 
jet events when higher order corrections are included. To start we consider first the cross section in 
Born approximation. Although not useful in practice, we may write the differential cross section in 
terms of jet-momenta by introducing delta-functions: 

dRniPa +Pb,Pl,-- ■,Pn)\M{pa,Pb,Pl,- ■ ■,Pn)f 
X6{Ji -Jiipi,.. .,Pn)) ... 6{Jn - JniPl, ■ ■ ■,Pn)), (2.3) 

where 5 denotes the center of mass energy squared, |A1p is the squared matrix element and dRn is 
the Lorentz invariant phase space measure 

n ,4 

dRniRpu .. .,Pn) = {2nf6{P - pd (2-4) 

i i=l '' ' 

where m,- denotes the mass of the /-th parton. The incoming particles with momenta pa and pb are 
assumed colorless. In case strongly interacting particles are considered in the initial state additional 
convolutions with the parton distribution functions need to be introduced. The momenta of the final 
state partons are given by pi,... ,p„. The functions Jiipi,. ■■,Pn) describe how the jet momenta are 
calculated from the parton momenta. Since in leading-order no recombination is possible, the jet 
momenta are identified with the parton momenta: 

Ji{pi,...,Pn) = Pi. (2.5) 

Obviously, it is than straightforward to evaluate the delta-functions and obtain the differential cross 
section in terms of the jet-momenta. Including next-to-leading order corrections, we need to con¬ 
sider the contribution from virtual corrections as well as real corrections. Ignoring for the moment 
the fact that both contributions are individually divergent due to soft and collinear singularities, we 
may apply the same argument as above to calculate the virtual corrections to the differential cross 
section. For the real corrections, however, the situation becomes more complicated. If we ask for 
precisely n jets in the final state, we need to integrate over the regions of the n + l parton phase 


dcr 1 r 

d^Ji...d^Jn 2s ^ 
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space in which n + \ partons are clustered to n jets. More precisely we need to evaluate integrals of 
the form 


jdR 

/l+l {Pa + Pb,Pi,- ■ .,Pn+l)\M{Pa,Pb,Pl,- ■ ■ ,Ph+i)P 
Xd(7i -Jiipi,...,p„+i))...6{J„-Jn{pi,... ,PH+l))0n-jet 


(2.6) 


where the functions J„ encode now how the n + \ partons are clustered to n jets and 0n-jet restricts 
the integration to the n jet region. (We note that the inclusion of the n + 1 jet region is straight¬ 
forward since no clustering occurs.) The functions /„ depend on the phase space region through 
the recombination procedure since in different phase space regions different partons are merged to 
form a jet. Evidently, the delta functions cannot be integrated numerically and an analytic approach 
is required. This is one facet of the problem we address in the next subsection. There is, however, 
a further problem: Using the standard recombination procedure, which is often simply the sum of 
the four momenta of the merged objects, we obtain in general massive jets even in case that we 
started with massless objects. As mentioned in the introduction it is thus a priori not clear how the 
contribution of the real corrections can be combined point-wise with the virtual corrections where 
the jets may have different masses. This is, however, required to define an evenf-weighf wifh NLO 
accuracy. In fhe nexf secfion we will show how fhe fwo issues are connecfed and can be addressed 
by a modificalion of fhe clustering prescripfion. In particular, we show how—by using a modified 
recombination procedure—fhe ‘real’ phase space can be faclorised info an n jef phase space and a 
remainder wifh fhe properfy fhaf fhe n jef phase space preserves fhe Born kinematics. As long as fhe 
fransfer function is nof approximated by a d-funcfion one could in principle relax fhe requiremenf 
fo map fhe unresolved regions of fhe real correcfions onto fhe Bom kinematics, since eq. (2.2) may 
be calculated for an arbifrary sef of parfonic variables used to describe virfual or real correcfions. 
However, using fhe aforemenfioned idenfificafion of fhe phase space for real and virtual correcfions 
allows to define poinf wise evenf weighfs in NLO accuracy. If is fhen sfraighfforward fo generalise 
fhe Mafrix Elemenf Mefhod fo NLO: The sef of y variables in eq. (2.2) are jusf reinterpreted as 
describing ‘fheory jefs’ as infroduced before. No furfher extension of eq. (2.2) is required. Note 
fhaf fhis approach is mean! fo be applied in fixed order: Parfon shower correcfions which parfially 
resum higher order correcfions would lead fo a double counfing if naively included in Ibis approach. 

2.2 A modified recombination prescription and phase space factorisation 

Using the resolution y,j to define the jet function F'j^ j (pi, .. .,Pn+\) we may write for the n-parton 
final state 



(2.7) 


The resolution depends on the final state objects: partons in case of the perturbative calculation and 
hadrons in case of the experimental analysis, ycut defines a preset value for the jet resolution. The 
momenta pi,...,p„ refer to the parfon momenta, while 7i,..., refer to the jet momenta. We have 


included them in the definition of the jet function since every jet algorithm includes in addition to 
the resolution also a prescription how to define the momenta of a jet in case it is formed by the 
merging of two finale state objects. In leading-order the jet momenta are identified with the parton 
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momenta, since no merging is possible. In NLO we also need the jet function for the case that 
n+\ partons form n jets. Since for soft and collinear configurations the jet function j needs 
to reproduce j to ensure the cancellation of soft and collinear divergencies, the jet function 
may be written as 

PTu...jSP^^ ■ ■ • - Z ©O'cut -y//Pl,.. • (2.8) 

i*j 

with 

&(ycut-yij(.pu---,Pn+i)) = 1 (2.9) 

for soft or collinear configurations. (For the moment we ignore initial state singularities. As we 
shall see the extension to include them as well is straightforward.) The step functions assure that a 
recombination of two partons into one jet occurs. For each possible combination i, j with yij < ycut 
the momenta 7, are obtained through the respective recombination procedure from the original mo¬ 
menta Pi- As mentioned before the mapping should respect momentum conservation and keep the 
recombined particle on the respective mass-shell in the sense defined above. This can be achieved 
by using fhe mapping infroduced in fhe Cafani-Seymour subfracfion mefhod. Depending on fhe 
unresolved partons and fhe chosen specfafor four differenf mappings are given in ref. [22, 23]. 
For each unresolved configurafion we may choose for example fhe combinafion wifh fhe smallesf 
momenfum fransfer fo fhe specfafor parfon. More general we may define funcfions 

(dk{p\,...,Pn+i),k = l,...,u-i-l (2.10) 

wifh fhe requiremenf fhaf 

2]0, = 1 (2.11) 

k 

fo selecf a specific mapping. For fhe numerical phase space infegrafion using Monfe Carlo mefhods 
if mighf be useful fo use smoofh funcfions 0^: insfead of sfep funcfions. 

Infroducing Yjk®k = I in eq- (2.8) we gef 

FXljSPU---,Pn+l)= Yj ®(ycut-yij{PU...,Pn+l))QkFl_j^iJu...,Jn) ( 2 . 12 ) 

For fhe jef cross secfion fhe confribufion from fhe real correcfions reads 

r/cr = ^ ^ @{ycul-yij{pu---,Pn+\))@kF\ j^^{J]_,...,Jn) 

X\Mn+\\^dRn+\{pi,. .. ,p„+l), (2.13) 

wifh fhe phase space measure as defined in eq. (2.4). In whaf follows we consider fhe funcfions 
@k as parf of fhe jef algorifhm. The role of fhe 0^ is fo selecf in each phase space region where 
parfons/jefs are merged fhe appropriafe clusfering. Since in each region @k selecfs a mapping 
(pi,... ,p«+i) —> (7i,...,/„) we may change fhe infegrafion variables accordingly using fhe respec¬ 
tive Cafani-Seymour parameferisafion of fhe phase space: 

<^cr = ^ @{ycvii-yij{P\,---,Pn+\))®kF\ j^{Ju---,Jn) 

X\M n+\ fdRn{Jl,...,Jn)dRij^k 


( 2 . 14 ) 


where dRij^k denotes the respective phase space measure introduced in ref. [22, 23]. Note that 
pi,.. .,Pn+\ appearing in the jet function and the matrix elements should be expressed in terms of 
the the momenta 7i,... and the integration variables used in dRij^k- Using the factorised phase 
space it is straightforward to calculate the contribution of the real corrections to the event weight 
in NLO accuracy: 

^ ^ ®(ycut-yij{Pl,---,Pn+l))Qk{Pl,---,Pn+l)\Mn+\\^dRij^k- (2.15) 

In the integration the momenta p\,...,pn+\ are determined from the jet momenta and 

the variables used in dRij^t- The inversion of the mapping (pi,...,p„+i) ^ (7i where <1) 

denotes the collection of variables used in the (unresolved) phase space measure dRtj^t is discussed 
in section 3. 

We have ignored so far that the phase space integration is in general divergent due to soft and 
collinear singularities. Since the Catani-Seymour phase space factorisation is valid in d dimension 
it is straightforward to regularise the divergencies within dimensional regularisation. Conceptually 
the easiest way to deal with the singularities is to apply a phase space slicing [25, 26]. In the 
numerical integration the integration over the unresolved parton is cut-off to avoid the collinear and 
soft configurations. In the singular regions, soft and collinear factorisation can be used to simplify 
the matrix elements such that the integration can be done analytically. The singularities obtained 
in this way are then combined with the virtual corrections. 

As far as the application of the Catani-Seymour subtraction method is concerned the situation 
is more involved: To allow the combination of the (integrated) subtraction term with the virtual 
corrections the jet algorithm or in general the observable is evaluated for the reduced kinematics in 
the Catani-Seymour formalism. The term which is added (and subtracted) thus reads: 

Y Xi P\,...jS^^'---'^>^^'^d,kdRn+\{p\,---,Pn+\) (2.16) 

where denote the dipoles defined in the Catani-Seymour subtraction method. Note that the 
mapping to obtain the jet momenta 7, from the parton momenta pi is encoded in the dipole. We 
are not free to chose the mapping in this case as this would result in a mismatch with the contribu¬ 
tion integrated analytically and combined with the virtual corrections. The contribution from the 
subtracted dipoles can thus not be combined point wise with the real corrections calculated using 
eq. (2.15). In ref. [20] a similar conclusion regarding the application of Catani-Seymour dipole 
subtraction within that method is drawn. 

So far we have assumed only final sfafe singularities. The above approach can be easily extended to 
initial state singularities. The jet function needs to be extended to cover also initial state singulari¬ 
ties. Using the different mappings as introduced in refs. [22, 23] is sufficient to handle all different 
cases. 
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3 Phase space parameterisation 


The parameterisation of the n + \ particle phase space in terms of an n particle phase space times 
an ‘unresolved’ phase space follows the phase space factorisation as given in the context of the 
Catani-Seymour subtraction method [22, 23]. As mentioned before, the mapping of the real phase 
space to the reduced kinematics defines the clustering prescription for the 3 —> 2 jet algorithm, 
generalising the 2 —> 1 clustering normally used. Since in the modified jet algorithm the resolution 
is not affected and the recombined jet momenta reproduce the naive soft and collinear limits^, the 
modified jet algorithm automatically fulfills infrared safely, factorisation of initial state collinear 
singularities, and momentum conservation while keeping the resulting jets on-shell. Furthermore, 
the phase space factorisation allows to span the respective real phase space associated with each 
point in the n-jet phase space in a straightforward manner. 

Each combination of unresolved partons i,j picked by the resolution of the jet algorithm (0(ycut - 

yijipi,.. .,pn+\)) = 1) and the spectator k selected through 0*: defines a specific mapping to cluster 

ij,k 

n+l partons to n jets (pi,...,p„+i)-> Depending on whether j and k are final or 

initial state particles (/ is always a final slate parton) there are four qualitatively different types of 
mappings which can be formulated for massless or massive particles [22, 23]. To apply the method 
outlined in section 2 we need to invert these mappings. For a given set of on-shell jet momenta 
(7i,.. .,J„) and a set of variables describing the unresolved phase space (O) we need the mapping 

(/i,...,7„,(l)) ^ {pi,...,pn+i) (3.1) 

to generate the n + \ parton phase space. In the following subsections we collect the required 
formulae. As mentioned in the introduction related formulae can be found in refs. [24, 27, 28]. 

3.1 Final state clustering with final state spectator 
Massless particles 

As described in ref. [22] the phase space of n -i-1 massless partons can be factorised in terms of a 
phase space of n massless momenta—which we identify with the jet momenta 7,—and the dipole 
phase space measure dRij^k related to the emission of an additional parton: 

dRn+l [pa+Pb,P\,- ..,Pn-2,Pj,Pk,Pi) = dRn {pa + PbJu-.., Jn-2, J j, Jk)dRij,k (3.2) 

with the n-particle phase space as defined in eq. (2.4). The incoming momenta are given by pa 
and ph. The dipole phase space measure as given in Ref. [22, (5.20)] reads in four space time 
dimensions 

dz dy{l -y)Q((p{2K-4>))Q{z{l -z))&(y(l -y)). (3.3) 

2{2n) 

^In the soft limit (pi —> 0) the kinematics used in the Catani-Seymour formalism reduces to {pi,.. .. ■,Pn+\ 1, while 

in the collinear limit (p; —> zp,Pj —> (1 -z)p) the set of momenta reduces to {pi,.. . • .,j^, •.. ,p«+i,p). 
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The set of variables O = {0,z, j) used to parameterise the phase spaee dRij^k is diseussed below. The 
phase spaee parameterisation eorresponds to the following elustering of n + 1 partons to n jets 

y 

Jj = Pi + Pj- - Pk, (3.4) 

i-y 

Jk = - Pk, (3.5) 

1 - 3 ; 

Jm = Pm, for mi^j,k (3.6) 

n 

whieh fulfill momentum eonservation (7, = P) and the respeetive on-shell eonditions (7? = 

i=l 

0 , i = To invert the mapping (pi,...,p„+i) ^ (7i,...,7„,<1>) we first observe that the 

momentum pt and the momenta pm {m + i, j, k) ean be obtained in terms of the momenta 7, and the 
variable y through 


Pm ~ Jm, (3.7) 

Pk = {f-y)Jk. (3.8) 

To determine the missing momenta p, and pj we use 

Pij = Pi + Pj= Jj + yJk, (3.9) 

^ij = {-^j + y-^kf = 2y(Jj ■ Jk). (3.10) 

The momenta p, and pj ean be easily expressed in the rest frame of pij rotated sueh that 7y points 
along the positive z-axis. Using 2(7y • pij) = stj the momentum 7' in this partieular frame is given 

by 

7'= ^(1,0,0,1). (3.11) 

The momentum J'j is obtained from the given 7y by a boost into the rest frame of p,y and two 
subsequent rotations to annihilate the x and y eomponent of Jf 

7 ;. - A'-^ict>,)A';i9y)A\pij)Jj, (3.12) 

with the Lorentz transformations A^J^cpx), Ay{6y), and A^(p,y) given in the appendix. We used 
the hat to denote the parity transform of a four veetor x: x = (x^,-^. The angles 9y and are 
determined through 


cos(0,) - ^ sm{9y) - ^ (3.13) 

^(7^2+ (7^2 ^(7^)2+ (7^2 

J(7p2 + (7p2 _jy 

cos{(px) = - 5 -, sin( 0 J ^ (3.14) 

l^yl \Jj\ 

In this frame the momenta p'. and p'. read 

^ J 

p'i - —^{^,2^Jzif-z)cos(|), 2^Jz(f-z)sm(f>, 2z-f), (3.15) 

p'j - Pi- (3.16) 
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where we have used the definition [22, (5.6)] 


_ l{pi • Jk) 

lUj-Jk)' 

The momenta Pi,Pj follow from p'-,p'; by inverting the Lorentz transformations: 

Pi ^ A^{pij)K'y{-0y)A’'J,-(t)^)p\, 

Pj = J^^{Pij)K{-Gy)Ki-^x)p'j- 

(The inverse of h^ipif) is given by PJ’ipij).) 


(3.17) 


(3.18) 

(3.19) 


Massive particles 


For massive partons i, j, k the phase space can again be factorised in terms of a phase space of n 
jets and the dipole phase space measure dRij^k related to the clustered parton [23] : 

1 (F, p 1 ,..., Pn -2 ,Pi,Pk,P j) = dRn Jn -2 , J j, dk) dRi j^k ■ (3.20) 


The n-jet phase space is again given by eq. (2.4), where some of the are non-zero now. The 
dipole phase space measure as taken from Ref. [23, (5.11)] reads in four dimensions 

q2 ( 1 -^ 2 _^ 2 _^ 2)2 




4{2nr 




-©(l -Pi-Pj-Pk) 


xd4) dz dy{In-4)))e{{z-z-){z+-z)) ©((y-y-)(y+ -y)) 
with the Kallen function defined by 


(3.21) 


A(x,y,z) - +y'^ +z^-2xy-2xz-2yz, 


(3.22) 


and 


nin 


l~^n ~ 


The integrafion boundaries are given by [23, (5.13)] 

2 piPj 


^ij ^ Jdj, Q = Pi + Pj + Pk = Jj + Jk- 


y- ^ 


2 ’ 


- 1 


^-pf-Pj-K 

2 pk{l -pk) 


1 2 2 2 ’ 


Z± - 


2P + {1 


■Pi -Pi 


-pI): 


2[pj+p^. + [l-pj-pypl)y] 

with the relative velocities between p,- -i- pj and p,- or pk [23, (5.14)] 


[l±Vij,iVij^k), 


{1-P-P-Ml)y + 2P 


V 


{i-p-p-fi 

)(l-.v)] 


( 

\-pI-p]-pI) 

(1-y) 


(3.23) 

(3.24) 

(3.25) 

(3.26) 

(3.27) 

(3.28) 
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The phase space parameterisation corresponds to the following clustering of n + 1 partons to n jets 
[23, (5.9)] 


Jk = 


4 

[Api,pf) 

r 

(i.fA) 


Pk + 


,. 2 ,. 


Q 

2 


Jj = Q-Jk, 


(3.29) 

(3.30) 

(3.31) 

n 

which again fulfill momentum conservation (X 7, = P) and and the on-shell conditions (7? = 


Jm = Pm, {tntj,k) 


i=\ 


m 


r., Jf - nq for I + j). To invert this clustering the momenta pk and {m + i,j,k) are calcu- 


lated first: 


Pk = 

Pm ~ Jm, 

where we have used 




A 


V 



+ [d - j) (1 --^2 _^2j ^ 2/|2] M, (3.32) 

(3.33) 


'JiPk ■Q) = Q^ + pI- {Pi + P jf = Q^ + ml- Sij, 


together with the definition 


3^ ^ 


2 {pi-pj) 


0^ -mJ - nr. - 


(3.34) 


(3.35) 


Similar to the massless case it is convenient to express p, and pj in the rest frame of pij = pi + pj 
rotated such that Q in the respective frame points along the positive z-axis. Using 2{Q ■ pij) = 
+ Sij - the momentum Q' in this particular frame is then given by 


Q’ = 




2 


a- 




2 


(3.36) 


Again Q' is obtained from Q through a boost to the rest frame of pij and subsequent rotations: 

Q' = K''Mx)K,{ey)K\pij)Q, (3.37) 

where the angles are similar to eq. (3.13) and eq. (3.14). In this frame the momenta p'. and p'. read 


Pi = 


Pj = 


+ P^ - pj],\Pi\{^i^0' cos^, sin0'sin0, cos0')j> 


{ 1 

f Sij 

\2V^' 


/ I 

! Sij 

\2V^' 

iG^ 


(3.38) 

(3.39) 


with |^;| - ^ipf)^-m] = ^A{sij,mlmj). Using 

Sik - {Pi + Pkf= Q^[z{l-y)(l-p]-p]-p^^+p]+p^^ 


(3.40) 
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one gets 


Q\\ - J)(l _ 2z)-^? +^ij] 

COS0' --^- - - , . — - - - 

Under the exchange yu^ yu^,z ^ 1 - z we have cos(d') —> -cos(d') as it should be. 

The momenta Pi,Pj follow again from p'j,p'j by 

Pi = A‘>{pij)A;i-ey)Ki-c!>,)p'., (3.42) 

Pj = ^''iPij)^yi-Oy)K(-^^)p'r (3-43) 


3.2 Final state clustering with initial state spectator 


Massless particles 


Using Ref. [22, (5.45)] the phase space can be factorised into the phase space of n massless particles 
and the dipole phase space dRij^a- 

dRn+\[pa + Pb,P\,---,Pn-l,Pj,P^ = dRn[xpa+Pb,J\,.--,Jn-l,Jj)dRij,a. (3.44) 

The dipole phase space measure in four dimensions is given by [22, (5.48)] 

dRij,a - -]z:^d(l)dzdxQ{(p{ln-4)))<d{z{l-zy)&{x{\-xy). (3.45) 

2{2ny 

Note that dRij^a includes an integration over x leading to a convolution of the measures given in 
eq. (3.44). The space parameterisation corresponds to the clustering of n + 1 partons to n jets 


Jj = Pi + Pj-{l-x)pa, (3.46) 

Jm = Pm (m + ij) (3.47) 

n 

which fulfills momentum conservation (X 7; = xpa + Pb) and the on-shell conditions (7? = 0, 

1=1 

I - l,...,n). Inverting this clustering allows to parameterise n-i-1 partons by means of the n jet 
momenta and three integration variables x,z,^ as follows 

Pm = Jm, (3.48) 


and because of 


Pij = Pi + Pj= Jj + (1 - x) Pa, (3.49) 

Sij ^ {jj + {l-x)Paf , (3.50) 

Pi and Pj can be calculated using the steps outlined in eq. (3.11) - eq. (3.19). 
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Massive particles 


Using ref. [23] the phase space of n + 1 massive partons can be expressed as a phase space of n 
particles convoluted with the dipole phase space dRij^a- 

dRn+l [pa + Pb,Pl,---,Pn-l,Pj, Pi) = dRn [xpa + Pb, J , Jn-lJ j) dRij,a- (3.51) 

The dipole phase space measure dRij^a reads [23, (5.48)]: 

^ d(f>dzdx&i(f>i27r-(p))&{(z-z-)iz+-z))&ixix+-x)). 


2{2n^) 

The integration boundaries are given by 

x+ = l+pjj-[pi+pjf, 


1-x + pI+p^-P^j± ^ 

l(^\-x + p 



2 ( 


\ 


dRiUa - dz dx @{(I){27 t- (l)))&({z-z-){z+ -z))&{x{x+ -x)). (3.52) 

(3.53) 

z± = - ^^^^(3.54) 

Pn - , / , mj = ^jJr (3.55) 

^^2Jj ■ Pa 

The phase space parameterisation corresponds to the clustering of n + 1 partons to n jets as in 
eq. (3.46) and eq. (3.47) but satisfying now the on-shell conditions = mf. and 7? = mj for I + 
j. To invert the mapping (pi,...,p„+i) ^ (/i,...,7„,0) (O = {x,z,^}), we start again in the rest 
frame of = pi + pj rotated such, that the momentum Jj points along the positive z-axis. Using 
{pij ■ J j) - j{sij - rn^j) the corresponding momenta 7'. in the rest frame of pij is given by 


with 


m„ 




The relation to 7y is again given by a sequence of one Lorentz boost and two rotations: 


J': - Ki<PM6y)A‘’ipii)Ji. 


the momenta p'. and p'. are given by 


Pi = 


Pi = 


t 2.2 

Sij-m.+mf 

2^j 

Sij-inf + nfi 


,Ip-I(sin0'cos0, sin0'sin0, cos0') 


2 

with|p;h = 


-Pi 


(3.56) 

(3.57) 

(3.58) 

(3.59) 


2 2 
O’ i ’ J 


) and 


mf-mf-(l-2z)Si! 
COS 0 ' =- ■' 


^A(sij,mlmj) 


(3.60) 


Pi and Pj follow from p'. and p'. by inverting the Lorentz transformations given in eq. (3.57). 
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3.3 Initial state clustering with final state spectator 


Massless particles 


The phase space of « +1 massless partons can be expressed as a phase space convolution of a phase 
space of n massless jets and the dipole phase space measure dRia^k for the emission of an additional 
massless parton from the initial state with a massless final state spectator. Most statements from 
section 3.2 can be carried over by the replacements a^k and y —> a (see Ref. [22]). However, since 
we now are dealing with clustering in the initial state, collinear singularities must be factorisable 
into the parton distribution functions. Because 


Pa-{Pi + Pk)-Pi-Pk 
Pa-{Pi+Pk) 






(3.61a) 


a jet function applying this clustering fulfills the condition for factorisability of initial state collinear 
singularities 


•^Pn-\,Pk,Pi', Pa,Pb) 


- j {P\,■.■,Pn-l,Pk^,ZPa,Pb)■ 

Pi^(^-z)Pa " 


(3.61b) 


Massive particles 

The phase space of n massive partons and one massless parton can be expressed as a phase space 
convolution of a phase space of n massive jets and the dipole phase space measure dRia^k for the 
emission of an additional massless parton from the initial state with a massive final slate spectator. 
All statements from section 3.2 can be carried over by the replacements k and j —> a, nii 0 
and niij —> mj, (see ref. [23]). The argument from eq. (3.61) also holds. 


3.4 Initial state clustering with initial state spectator 

In case of initial state clustering with an initial state spectator the phase space can again be written 
as a convolution [22, (5.149)]: 


dRn+\ {pa+Pb,Pl,---,Pn,Pi)^dRn {xpa + Pb, J I, ■ ■ ■, J n) dRia,b, (3.62) 

with the n-particle phase space given in eq. (2.4). The dipole phase space measure dRia^t reads [22, 
(5.151)]: 


dRiab - ^^ dcf) dv dx 0 (0(27r - 0)) 0(t;) 0(1 -— )0(x(l -x)). (3.63) 

lilnY V l--r/ 

The phase space parameterisation corresponds to the following clustering of n +1 (massless/massive) 
partons to n (massless/massive) jets: 

Jm= ^ia,bPm, m=l,...,n (3.64) 


with 


Pa+Pb-Pi, K= XPa + Pb, 


(3.65) 
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and the Lorentz boost transforming K into K given by [22, (5.144)]: 


{J^ia,bfy = g^v ■ 


l(K + Kf(K + K\ ^ 2K‘'K, 


{K + Kf 




(3.66) 


The inverse boost is obtained by exchanging K and K. All outgoing momenta pi are transformed 

n 

to balance the transverse momentum. Momentum conservation (!]// = xpa + Pb) and on-shell 

i=\ 

conditions - mj, Z =!,...,«) are not affected by the boost. Inverting this clustering allows to 
parameterise n -i-1 partons by means of the n jet momenta and three integration variables x,v,(p as 
follows. The momenta p,„{m- I,...,n) are obtained by inverting the boost: 



II 

> 

(3.67) 

Using the definition for v 

Sia ^ VSab, 

(3.68) 

together with 

Sib ^(I-X-V) Sab 

(3.69) 

which leads to 

Sia + Sib x)Sab 

(3.70) 


it is straightforward to express the momentum pi in the rest frame of pa + pb rotated such that pa 
points along the z-axis. In this particular frame the momentum p, is given by 


p'i ^(1--^) 

^^^(l,sing-cos0, sin0; sin0, cos0-). 

(3.71) 

Using Pa in this particular frame 

^(1.0,0,1), 

(3.72) 

the angle B' can be read off 

, 2v 

COS0,- = 1-. 

' l-x 

(3.73) 


The momenta p, is obtained according to 


P, = A;(-0,)A;(-0Jp; (3.74) 

as in the previous cases with Jj —> pa- Note that no massive case as in the previous sections needs 
to be studied since the two incoming partons and the collinear parton are always assumed to be 
massless. 
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4 Consistency checks 


To validate the approach, we apply the procedure outlined in the previous sections to two processes. 
As an example, where one has to deal with initial state singularities, we study Drell-Yan production 
in hadronic collisions. More precisely, we calculate Drell-Yan production at the LHC running at a 
center-of-mass energy of = 13 TeV. We apply phase space cuts similar to what is used in the 
LHC experiments. For simplicity, we veto any additional jet in the final state since we are only 
interested in the case where recombination occurs in the real corrections. The inclusion of the 
contribution due to an additional jet is straightforward since no recombination occurs. For the final 
sfafe elecfrons, we require fhe invarianf mass nige of fhe elecfron pair in the region defined by 116 
GeV< ruee < 3 TeV. Furfhermore, we demand a minimum fransverse momenfum p> 25 GeV for 
fhe elecfron and resfricf fhe rapidify of fhe elecfron fo \'qe\ < 2.5. Unresolved inifial sfafe radiation 
is clusfered wifh fhe beam according fo secfion 3.4. 

As a second example, where one has fo deal wifh final sfafe radiafion, we analyse fop-quark pair 
producfion in e^e~ annihilafion. Similar fo fhe Drell-Yan case, we vefo again fhe emission of an 
addifional jef. For fhe fop-quark mass we use m, = 174 GeV. For fhe cenler-of-mass energy we 
choose ^fs = 500 GeV relevanf for a fulure linear collider. We do nol include fhe decay of fhe lop 
quarks, inslead, we Ireal Ihem as fagged lop-jels. These jefs are obfained wifh a kt-jel algorilhm 
wifh fhe resolufion criferia defined by 

yii = ^ -, (4-1) 

and the resolution ycut set to ycut = 0.1. For the recombination of unresolved particles the modified 
3 —> 2 clustering prescription according fo secfion 3.1 is used. 

Allhough very simple, Ihese Iwo examples cover essentially all relevanf cases. Furfhermore, com- 
pacl analylic resulls are available for fhe higher order corrections and if is slraighlforward fo apply 
fhe ideas ouflined in fhis arlicle. For delails on fhe NLO calculalions using phase space slicing we 
refer fo refs. [36, 37]. 

Exclusively demanding n jefs in fhe final sfafe allows fo define a differenlial n-jef even! weighl al 
NLO accuracy: 

j^NLO ^ ^ j^V ^ 

d^J\ ... d'^Jn d^j\ ... d^Jn d‘^Ji ... d'^Jn d^Ji .. . d^Jn 

We use fhe superscripls B, V and R fo indicale fhe confribulions from fhe Born malrix elemenls, 
fhe virlual correcfions and the real corrections. In case of the real corrections a regularisation 
of the soft- and collinear singularities using the phase space slicing method is understood. The 
‘unresolved’ contribution is included in the virtual corrections and cancels the respective soft and 
collinear singularities. We note, that the real corrections are calculated using eq. (2.15) which 
means that for each phase space point (7i,...,7„) an additional three dimensional integration is 
required to obtain dcr^. To check the approach and the numerical implementation we use eq. (4.2) 
integrated over the phase space to calculate the total cross section. The results can be compared with 
the ones obtained by a standard parton level Monte Carlo. We have checked that the results using 
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eq. (2.15) are in perfect agreement with the results of a conventional parton MC. One may argue 
that the comparison of the total cross section is not very sensitive to the details of the calculation 
and inconsistencies in specific phase space regions could escape detection. In fact eq. (4.2) can 
also be used to calculate arbitrary distributions: 





U Jffi 


r/crNLO 

d^Jl ...d^Jn 


d(0{Ji 


Jn)-0{Ji 



(4.3) 


Again these contributions can be compared with the outcome of a parton level Monte Carlo. This 
comparison allows a detailed check of the entire phase space. We stress that in the parton level 
Monte Carlo the same modified jef algorithm (3 —> 2 clustering!) has to be used. In figure 1 
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Figure 1. Differential distributions for Drell-Yan production calculated using a conventional parton level 
MC compared with a calculation using the factorised jet phase space as described in section 3.4. 


we collect various distributions calculated for Drell-Yan production. In particular, we show the 
angular and the energy distribution of the scattered electron/positron. In addition, the invariant 
mass distribution and the rapidity distribution of the system are given. The blue solid lines 
show the results obtained with a conventional parton level Monte Carlo. The red dashed lines 
show the results using the factorised jet phase space as illustrated in the previous sections. In the 
lower part of the plots we show the discrepancy between the two approaches in terms of standard 
deviations, where the uncertainty of each approach is due to the limited statistics of the Monte Carlo 
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integrations. For the parton distribution functions we use the CTlOnlo pdf set [38]. The center-of- 
mass energy is set to 13 TeV and we applied the aforementioned cuts. For all four distributions, we 
find perfect agreement between the two approaches. In most cases the discrepancy is less than one 
standard deviation. In figure 2 a similar analysis is shown for top-quark pair production in e^e'~ 
annihilation. In particular, we show distributions with respect to the cosine of the azimuthal angle 
of the outgoing top quark, the polar angle distribution, the transverse momentum distribution and 
the rapidity distribution. Again the blue solid curves show the results of a conventional parton level 
Monte Carlo while the red dashed curves give the results using the factorised jet phase space. Note 
that in both cases the modified 3^2 clustering is employed. Again we find perfect agreement 
between the two approaches. 
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Figure 2. Differential distributions for top-quark pair production in e^e~ annihilation calculated using a 
conventional parton level MC compared with a calculation using the factorised jet phase space as described 
in section 3.1. 


4.1 Impact of NLO corrections — k-factors 

It has been pointed out in ref. [18] that restricting the NLO analysis to the Bom level kinematics 
may lead to rather moderate k-factors in general. In figure 3 we show the respective k-factors for 
the previously studied differential distributions. In case of the angular distribution of the outgoing 
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Figure 3. Impact of NLO corrections on differential distributions for Drell-Yan production at a hadron 
collider. 


electron, NLO corrections at the level of only a few percent are observed. For the rapidity distri¬ 
butions they are slightly larger but still small in absolute size. For the energy distribution of the 
electron and the invariant mass distribution of the lepton pair the corrections seem to be larger. 
However, the ^-factor suffers from statistical uncertainties and shows large fluctuations. In regions 
where the statistical fluctuations are small we find again a moderate A:-factor. In figure 4 the A:-factor 
for top-quark pair production in e^e~ annihilation is shown. In all cases we find NLO corrections 
of a few per cent only and thus k-factors very close to one. We thus extend the observations of 
ref. [18] also to final state radiation. 

4.2 Impact of modified clustering / jet algorithms 

All the previously shown differential distributions have been obtained using the modified jet algo¬ 
rithm: in the conventional parton level Monte Carlo as well as in the alternative approach using 
a factorised jet phase space. As pointed out in the introduction we consider the modification of 
the clustering as part of the intrinsic ambiguities of jet algorithms. As a consequence we do not 
expect a large effect of the clustering. If in contrast a large effect is observed one should question 
the definition of the observable since it shows a large sensitivity on an aspect of the jet algorithm 
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Figure 4. Impact of NLO corrections on differential distributions for top-quark pair production in e'^'e' 
annihilation. 


which is not well defined. Since it is difficult to make general statements about the size of possible 
effects, one should investigate the impact of the new clustering on a case by case study to make 
sure that no large deviations are observed. In figure 5 the conventional 2 —> 1 clustering which does 
not respect the on-shell condition of the clustered objects is compared for Drell-Yan production 
with the 3 —> 2 clustering advocated here. Note that both results have been obtained using a con¬ 
ventional parton-level MC. The blue solid curves show the result using the 2 —> 1 clustering, while 
the red dashed lines give the results for the 3 —> 2 clustering. Since for Drell-Yan production the 
clustering never includes the outgoing electron we do not expect a major effect. Indeed figure 5 
shows essentially no difference within the statistical uncertainty. A minor effect is visible in the 
angular distribution and in the rapidity distribution. This can be related to the initial state cluster¬ 
ing which may introduce an additional boost orthogonal to the beam axis which can influence the 
angular distributions. In figure 6 the corresponding result is shown for top-quark pair production. 
For the angular distribution the two different algorithms give the same result within the statistical 
uncertainties. For the transverse momentum distribution a large effect is visible at large transverse 
momentum. This is not surprising since at phase space boundaries we expect to become sensitive 
to the details of the clustering. Below 160 GeV we observe that the 3^2 clustering leads to dis¬ 
tributions which are between two and five per cent below the traditional 2 —> 1 combination. In 
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Figure 5. Impact of 3 —» 2 clustering with respect to 2 —» 1 clustering on differential distributions for Drell- 
Yan production at a hadron collider. 


general the 2 —> 1 clustering leads to an increase of the mass of the clustered object which could 
be responsible for the observed pattern. To analyse this effect we show in figure 7 the distribution 
of the mass of the top-quark jet using the conventional 2 —> 1 clustering. As one can see most of 
the events have a jet-mass close to the nominal top-quark mass. However, there are also events 
with jet masses up to 280 GeV. Note that using the modified clustering the jet mass is fixed to the 
top-quark mass. In particular at phase space boundaries, the difference in the jet mass may result in 
distortions of distributions which are sensitive to mass effects. This is precisely what we observe in 
the lower plots of figure 6 where one can see that indeed the largest effects arise at the phase space 
boundary. This is not surprising since minor changes in the mass of the clustered objects become 
important. Since the distributions are normalised this effect introduces also a modification in the 
distribution away from the phase space boundary. Using cuts to avoid the phase space boundaries 
should thus result in smaller differences between the two different clustering prescriptions. This is 
illustrated in figure 8 where we show the same distributions but now using additional cuts. Again 
the blue solid line shows the conventional 2^1 clustering while the red dashed line shows the 
alternative 3 —> 2 clustering. Indeed we find that the difference becomes smaller and is of the order 
of 1 % only, which might be seen as an intrinsic uncertainty. To close this section we stress that 
minor differences between the two clusterings are not per se problematic as long as everything is 
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Figure 6. Impact of 3 ^ 2 clustering with respect to 2 ^ 1 clustering on differential distributions for top- 
quark pair production in e^e~ annihilation. 
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Figure 7. Distribution of the mass of the top-quark jet using the conventional 2 —> 1 recombination. 


done consistently and the same clustering is used in the experimental analysis. 


- 24 - 














































































Figure 8. Same as figure 6 but with additional cuts to avoid the phase space boundaries. 


5 Application 

In this section we apply the MEM to top-quark pair production in e"^e~ annihilation. First we use 
the aforementioned procedure to generate unweighted events with NLO accuracy. These events are 
then analysed using the MEM in EO and NEO accuracy. In particular, we illustrate the extraction 
of the top-quark mass from the event sample. 


5.1 Generating unweighted jet events 


For a chosen jet algorithm with a preset value of the resolution ycut it is now straightforward to 
generate unweighted jet events using an ‘acceptance-rejection’ algorithm. The respective NEO jet 
weight is given by 


pU\ 


Jn) = 


d^J\ ... d^Jn 


(5.1) 


where the right hand side is evaluated according to eq. (4.2). The acceptance-rejection method 
requires an upper boundary p^ax of the weight This can be obtained for example 

within a phase space integration. An n-jet candidate event is then constructed using (in - 4) random 
numbers. As a measure for the probability the weight introduced in eq. (5.1) is calculated for the 
candidate event. Note that a three dimensional integration must be performed to do so. Generating 
an additional uniformly distributed random number ru between 0 and p^ax the candidate event is 
accepted if ru is below the aforementioned weight. 

In principle it is also possible to generate unweighted NEO n-jet events (7i,...,/„) together with 
n+ 1-jet events (/',...,7'^^) from n + \ partons (pi,...,p„+i) by augmenting the definition of p: 


piPi 


Pn+l) = 


dcr 


.NLO 




d^Jl.-.d^Jn 

Jo-NLO 


«+l 


(5.2) 
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The let functions F , and F ,, decide whether n or n + 1 lets are resolved and how the 
momenta pi are clustered into the jets. The n + 1-jet events (/', are obtained by the identi¬ 

fication 

J'^Pi. (5.3) 

The n-jet events (7i,7„) follow from clustering by the 3 —> 2 jet algorithm 

{pi,...,Pn+l) iJl{pi,...,Pn+\),...,JniPU-,Pn+l))- (5.4) 


The main difference wifh respecf fo fhe previously described even! generation is, fhaf now n + 
1 parfon momenfa are generated using (3(n-i- l)-4) random numbers. While fhis mefhod will 
generafed n-jef evenfs wifh NLO accuracy we sfress fhaf fhe generafed n + 1-jef evenfs have only 
LO accuracy. 


To validate fhe generation of unweighfed evenfs, we reproduce fhe differential disfribufions calcu¬ 
lated in section 4. In fofal we generafed 73128 evenfs wifh NLO accuracy. As in section 4 we veto 
fhe emission of an addifional jef. In figure 9 we show fhe comparison wifh disfribufions calculafed 
using fhe conventional parfon level Monte Carlo. The blue solid lines represenf fhe resulfs from 
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Figure 9. Validation of the generation of unweighted NLO top-quark pair events (no additional jet) 


the parton-level Monte Carlo while the red dashed lines show the distributions calculated from the 
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unweighted jet events generated as described above. Below we show the difference between the 
two distribution in units of one standard deviation. As one can see we find perfecf agreemenf wifhin 
fhe sfafisfical uncerfainfies. 

5.2 Matrix Element Method 

The possibility to generate unweighted jet events with NLO accuracy together with the possibility 
to assign NLO event weights to them, allows to perform a validation of the Matrix Element Method 
at NLO using the generated events as input to a toy experiment. As a concrete example we illustrate 
the extraction of the top-quark mass in e^e~ —> tt employing the MEM at NLO. We note that this 
study may be relevant for the top-quark mass measurements at a future linear collider. 

As mentioned before we ignore for simplicity the top-quark decay and assume that top-(anti)quark 
jets are observed. An event is than defined by the energies and angles of the respective jets 
{Et, cos Of, Ej, cos 9f, (pi)- This fixes the jet momenta depending on the top quark mass m, 
as 


Jt = {Et, |pf|cos(^fSin0f, |p;|sin0(Sin0f, \pt\cos9t), 

Ji = {Ej, \pi\cos(pism9i, Ipr|sin0fsin0f, \pi\cos9i), (5.5) 


with \ptj\ = ^E^--m^. Exclusively demanding a top-quark pair without an additional jet fixes the 
energies Et = Ej = V^/2. A 2-jet NLO event weight for x- {cos 9t, cpt, cos 9i, (pi) can be obtained 
according to 


dcrNLO 

dx 


d cos 9t d(pt d cos 9i d(pt 
pt dcT^^^ 


32n-2 d^Ji d^J2 


COS01 =COS0/, COS02=COS0f, (pl=4>t, 


with Pt = "y 1 —A sample of N unweighted 2-jet NLO events 

{4 = {cos9t, (p't, cos9\, (p\), i = 1,---,a} (5.6) 

is generated for some ‘true” top-quark mass = 174 GeV. The NLO likelihood for the 
sample can be constructed from the differential 2-jet cross section as follows 


N 


£^^^{mt) = Y]£^^^{xi\mt) = 

P( 


i=l 


1 


q-NLU(^^) 


^ dcr^^^{mt) 


i=l 


N N 


32n^ cr^^^{mt)) j i \d^Jt d^Ji 


n 


dcr 


.NLO 


dxi 

{mt) 


(5.7) 


where the dependence on rut is shown explicitly and Jt and Jj follow from x, according to eq. (5.5). 
Note that the jet momenta when evaluated for the event xi depend on the mass nit. The negative 
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logarithm of the likelihood (or ”Log-likelihood”) therefore reads 


N 


N 


^ A^log(cr^^°(mf))-^log 


(=1 


i=l 


dxj 


- A^log 


Jt 


N 


Yj^og 


i=l 


da^^^init) 
d^Jt d^Jt 


(5.8) 


event i 


Maximising (minimising) this likelihood (Log-likelihood) with respect to m, yields an estimator mt 

nil 

-log£^^° (nit) = inf(-logi:’^^°(mf)). (5.9) 

nil ^ ' 


The lefthand plot of figure 10 illustrates the Log-likelihood for 73128 events generated with NLO 
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Figure 10. NLO and Born Log-likelihood for 73128 NLO events between m, = 100 GeV and m, = 240 GeV 
(”trae” value mj™ = 174 GeV) and zoomed in between m, = 170 GeV and rtit - 180 GeV to extract nit and 
Am, by a fit 


accuracy as function of the top-quark mass mt in the range between 100 GeV and 240 GeV. The 
solid blue line shows the result evaluating the likelihood using LO predictions, the dashed red 
line shows the result obtained using NLO predictions. As mentioned before a top-quark mass 
^true _ GeV has been used to generate the events. The righthand plot of figure 10 shows Ihe 
area around Ihe minimum belween m, = 170 GeV and nit - 180 GeV in order lo exlracl nit and bSht 
by a parabola hi (see ref. [39]). As one can see from figure 10 exlracling Ihe lop-quark mass wilh 
a likelihood based on Ihe Bom approximation yields an estimator which shows a significanl 
deviation from Ihe inpul value m]™® hence can nol be regarded as an unbiased estimator. More 
precisely we find 

(178.7 ±1.2) GeV (5.10) 

which is 4 cr way from Ihe Irue lop-quark mass used in Ihe evenl generation. Note lhal in LO 
Ihe resulls are independenl of Us- It is thus not possible to attribute a theoretical uncertainty by 
simply varying the renormalisation scale. On the other hand extracting the top-quark mass with 
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the likelihood based on NLO predictions results in an estimator which is consistent with the 
input value nr*™® within the uncertainty. Using NLO accuracy we find 

(5.11) 


in perfect agreement with nr*™®. We stress that in both cases we use the same unweighted events 
generated with NLO accuracy. In principle the discrepancy between and nr^^® is not surpris¬ 
ing since we used in LO predictions to analyse events generated with NLO accuracy. What 
is however remarkable is the large size of the effect. As has been illustrated in section 4 the NLO 
corrections are usually small for most of the distributions. Nevertheless we observe a large effect 
using LO or NLO predictions within the Matrix Element Method. From the above results we may 
conclude that the Born matrix element evaluated for m, = 178 GeV gives a better approximation 
of the NLO corrections evaluated for m, = 174 GeV than the Born approximation evaluated for 
174 GeV. To investigate this point further we show in figure 11 fhe comparison of the two predic¬ 
tions. Obviously the NLO corrections cannot be completely absorbed by changing the mass in the 



cos 9t 


Figure 11. Comparison of the NLO predictions (red dashed line) evaluated for nit - 174 GeV with the Born 
approximation (blue solid line) evaluated for m, = 178 GeV. In the lower plot the ratio of the two is shown 
(black line). For comparison the k-factor for nit - 174 GeV is shown in gray (cf. figure 4). 


LO predictions. Comparing however the black with the gray line in figure 11 we find thaf indeed 
rrit = 178 GeV gives a slightly better description of the NLO result. In the range -1 < cos Ot < 0.75 
the difference is below 1% and fhe maximal deviation at cosOt = 1 is 4%. The difference is below 
1% in the range -1 < cos6? < 0.38 and the maximal deviation at cosOt = 1 is 6% when nit = 174 
GeV is used in the Bom approximation. 

In view of a future linear collider we stress that the renormalisation scheme is well defined in the 
above procedure. Applying the above procedure to realistic data, the top-quark mass within the 
pole mass scheme would be determined. 
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Figure 12. Distribution of the estimator rht around the true value = 174 GeV and dependence of the 
Log-likelihood on the number of events. At the bottom of the right plot the respective estimator m, + Am, is 
shown. 


In figure 12 we show the consistency of the approach. In particular, we illustrate that provides 
indeed an unbiased maximum likelihood estimator (ref. [40]). The lefthand plot of figure 12 shows 
the distribution of the estimator if we interpret our event sample with 73128 events as 203 
independent toy experiments with 360 events each. The dashed line shows a gaussian fitted to the 
data. The righthand plot illustrates how increasing the number of events results in a tightening of 
the dip in the Log-likelihood around the true value of m, and the approximate scaling of the error 
of the estimator Ain, oc N ~2 (see bottom of righthand plot). 

As a final remark we comment on the impact of the modified jet algorithm. Top-quark pair produc¬ 
tion in e'^e~ annihilation is highly constrained through momentum conservation and the underlying 
symmetries of the interaction. Most of the sensitivity to the top-quark mass stems essentially from 
information contained already in in the cos(0,) distribution. On the other hand for this distribution 
the two clustering algorithms give the same result at the permille level as we have shown in figure 6. 
As a consequence we do not expect major differences in case the conventional clustering would be 
used in the experimental analysis. 


6 Conclusion 

In this article we have shown how to calculate event weights for jet events at NLO accuracy. The 
ability to define event weights at NLO is a necessary prerequisite to extend the Matrix Element 
Method beyond the Born approximation. The basic ingredient of the method presented here is a 
modification of the clustering prescription used in jet algorithms. Instead of using the conventional 
2 —> 1 clustering, where the momentum of the clustered object is just the sum of the two initial 
jet candidates, we use a recombination inspired by the phase space mapping used in the Catani- 
Seymour subtraction method. This leads naturally to a factorisation of the phase space for the real 
corrections into resolved and unresolved contributions. Furthermore, the factorisation allows to in- 
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tegrate numerically the contribution of the unresolved configurations after an appropriate regulator 
to handle the mass and soft singularities has been chosen. Similar ideas have been investigated in 
a different context already in refs. [24, 27, 28]. The major difference is that we consider the new 
clustering not only as a technical trick but as an essential part of a modified jet algorithm. Using 
this modified clustering no artifical jet mass is generated and it is straightforward to map the events 
obtained onto the bom kinematics. As validation of the proposed method, we have successfully 
reproduced differential distributions at NLO accuracy in Drell-Yan production and top-quark pair 
production in e^e~ annihilation. Although simple, these two examples cover essentially all relevant 
cases. We have also investigated the impact of the modified jet algorithm. At phase space bound¬ 
aries the effects can be large. Additional cuts can be used to reduce the impact. The remaining 
effect may be considered as an intrinsic uncertainty inherent to jet algorithms. For the examples 
studied here the effect is reduced to the per cent level, after applying cuts. We stress however that 
in hadronic collisions the situation could be different and one needs to investigate the impact of the 
new clustering on a case by case study. As a further application we have studied as a toy exam¬ 
ple the Matrix Element Method at NLO applied to top-quark pair production in e"^e~ annihilation. 
More precisely, we investigated the determination of the top-quark mass. This study is relevant 
for a possible future Linear Collider. Applying the Matrix Element Method to events generated 
with NLO accuracy we observe that the MEM in LO fails to reproduce the input value. While the 
NLO analysis correctly reproduces the input with an uncertainty of about 1 GeV for about 70000 
simulated events, the LO analysis leads to a value off by 4 GeV. These findings should be taken into 
account, when top-quark mass measurements at the Tevatron using the MEM are discussed. Let us 
end with a final remark concerning parton shower corrections. As mentioned in section 2 the naive 
inclusion of corrections due to the parton shower would lead to a double counting. Lurther studies 
are required to extend the method presented here in this direction. 

Note added; 

While we were in the process of writing this article ref. [41] appeared, where an extension of the 
jet algorithm in e^e~ annihilation to massless quarks similar to what is discussed here has been 
presented. 
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A Explicit form of the Lorentz transformations 


For a given four vector X with 0 the rotational free boost from the rest frame of X to the 
system where X takes the form as given reads when applied to the four momentum y: 


A\X)y = 


X^ 


0 , (X-f) ,, 
=3^ +^^’3^ + 


^/x2 


(X-f) 


+ V^) VF. 




(A.1) 


Defining X = (X°, -X) the boost from the frame in which X is given to the rest frame is given by 
A^{X). In fact, eq. (A.l) is a special case of the more general boost given in eq. (3.66): 

a\X) = Aia,b, for K = i ^/^A K = (A°,X). (A.2) 


The Lorentz transformations for rotations around the x and the y axis are given by 


Ki4>) = 


= 


MO 0 0 

0 10 0 

0 0 cos(0) sin(0) 
0 0 - sin(^) cos(0) 

10 0 0 

0 cos(d>) 0 -sin(0) 

0 0 10 

0 sin(0) 0 cos(d>) 


(A.3) 


(A.4) 


For the product we have 


A;(0)A^(0) 


MO 0 O' 

0 COS0 sin^sin^ -sin^cosi;^ 
0 0 cos<p sin0 

,0 sin0 -cos6sin^ cos0cos0 , 


(A.5) 
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